# Replication code for Newell, Pizer, Prest, "A Discounting Rule for the Social Cost of Carbon"
# Code Author: Brian Prest, RFF
# Run on R version 4.0.2
rm(list=ls()); gc()
library(data.table)
library(lubridate)

# Useful functions
"%&%" <- function(x,y) paste0(x,y)
rep.row<-function(x,n) matrix(rep(x,each=n), nrow=n)
rep.col<-function(x,n) matrix(rep(x,each=n), ncol=n, byrow=TRUE)
getsize = function(object, scaling=2^20) {
  # Returns a vector of object sizes, including the total, in MB (2^20)
  obs = sort( sapply(object,function(x) object.size(get(x)) ))/scaling
  tot = sum(obs)
  return(round(c(obs, total=tot), 1))
}

root='/path/to/replication/dir/'
code.dir = root%&%'/src/'
data.dir = root%&%'/Data/'
figure.dir = root%&%'/Figures/'

# Set palette
palette(c(rgb(4, 39, 60, maxColorValue=255), # RFF black
          rgb(136,196,244, maxColorValue=255), # RFF blue
          rgb(255, 102, 99, maxColorValue=255), # RFF coral
          rgb(80,177,97, maxColorValue=255), # RFF green
          rgb(116,100,94, maxColorValue=255), # RFF brown
          rgb(118, 94, 165, maxColorValue=255), # RFF purple
          rgb(244, 162, 95, maxColorValue=255), # RFF orange
          rgb(235, 211, 103, maxColorValue=255), # RFF yellow
          rgb(224, 228, 231, maxColorValue=255))) # RFF light gray

# Read in aggregated global MSW data
global_growth_percap = read.csv(file=data.dir%&%'/aggregated_MSW_growth_data.csv')
rownames(global_growth_percap) = global_growth_percap$X # make row names the year
global_growth_percap = global_growth_percap[,-1] # drop year column
global_growth_percap = as.matrix(global_growth_percap) # convert to matrix
colnames(global_growth_percap) = 1:ncol(global_growth_percap) # simplify column names

# Make figure 1 - Growth distribution over time
source(code.dir%&%'/plot_growth_distribution.R')

# Generate BR Term structures at each starting rate
source(code.dir%&%'/BR_term_structures.R')

# Read in CE rates from the literature
source(code.dir%&%'/read_in_CE_rates.R')

# Calibrate rho and eta, and make figures 2 and A.2
source(code.dir%&%'/calibrate_rho_eta.R')

# Make figure 3 - stylized discounting paths
source(code.dir%&%'/make_stylized_figure.R')

# Run SCC calculations and make figures 4 and 5
source(code.dir%&%'/run_scc_calcs.R')

# Calibrate rho and eta for extra BR Term structures for more starting rates, make figure A.1
source(code.dir%&%'/extra_BR_term_structures.R')

